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Abstract 


The  Scan-MUSIC  (Multiple  Signal  Classification),  or  SMUSIC,  algorithm  was  developed  by  the 
Millimeter  Wave  Branch,  Sensors  and  Electron  Devices  Directorate,  Army  Research  Laboratory 
(ARE).  The  algorithm  improves  angular  resolution  for  target  detection  with  the  use  of  a  single 
rotatable  sensor  scanning  in  an  angular  region  of  interest.  This  algorithm  has  been  adapted  and 
extended  from  the  MUSIC  algorithm  that  has  been  used  for  a  linear  sensor  array.  Experimental 
and  computer  simulation  results  have  resolved  two  closely  spaced  point  targets  that  exhibited 
constructive  interference,  but  not  for  the  targets  that  exhibited  destructive  interference.  The 
interferences  occurred  due  to  the  relative  phases  of  their  returned  signals.  Therefore,  there  were 
some  limitations  of  the  algorithm  for  the  point  targets.  In  this  paper,  I  extend  the  application  of 
the  SMUSIC  algorithm  to  a  problem  of  complex  scatterer-type  targets,  which  is  more  useful  and 
of  greater  practical  interest.  The  SMUSIC  simulation  results  show  that  the  algorithm  can  resolve 
centroids  of  the  complex  scatterer-type  targets  in  which  they  are  close  together  within  the 
beam  width  of  a  K„  radar  antenna  with  frequency  diversification.  To  make  it  easier  to  see  the 
results  of  the  simulation  and  to  demonstrate  the  SMUSIC  algorithm  for  the  complex 
scatterer-type  targets,  I  developed  a  graphical  interface  tool  to  run  and  analyze  the  SMUSIC 
algorithm.  The  tool,  which  is  written  in  the  MATLAB  language,  allows  an  end  user  to  analyze 
the  performance  of  the  SMUSIC  algorithm  by  entering  a  combination  of  parameters  necessary 
for  the  simulation. 
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1.  Introduction 


Scan-MUSIC  (Multiple  Signal  Classification),  or  SMUSIC,  algorithm  was 
developed  by  the  Millimeter  Wave  Branch,  U.S.  Army  Research  Laboratory 
(ARL),  to  improve  the  resolution  of  two  closely  spaced  point  targets  using  a 
scanning  radar  antenna.  These  targets  were  within  a  beamwidth  of  the  an¬ 
tenna.  There  were  some  limitations  of  the  algorithm  for  two  point  targets 
based  on  the  experimental  and  computer  simulation  results  shown  in  other 
documented  research  [1-4].  The  algorithm  resolved  the  targets  that  exhib¬ 
ited  constructive  interference;  it  could  not  resolve  the  targets  that  exhibited 
destructive  interference.  In  this  paper,  I  extend  the  application  of  the  algo¬ 
rithm  to  a  complex  scatterer-type  target  problem.  The  results  show  that  the 
algorithm  resolves  the  centroids  of  complex  targets  that  are  within  the 
beamwidth  of  a  K^,  radar  antenna  with  frequency  diversification  regardless 
of  interference.  These  targets  consist  of  many  point  scatterers;  each  scat- 
terer  has  its  amplitude  and  its  own  phase  associated  with  its  returned  sig¬ 
nal.  I  have  developed  a  graphical  interface  tool  to  run  and  analyze  the  algo¬ 
rithm  for  the  complex  scatterer-type  targets.  This  tool,  written  in  MATLAB”^, 
will  enable  us  to  analyze  the  performance  of  the  MUSIC  algorithm  in  some 
practical  scenarios  and  will  help  to  improve  future  U.S.  Army  radar  systems. 

This  paper  is  divided  into  four  remaining  sections.  Section  2  describes  the 
extension  of  the  MUSIC  algorithm  to  the  SMUSIC  algorithm.  Section  3  de¬ 
scribes  the  interface  program.  Section  4  shows  how  to  use  the  interface  tool 
and  the  SMUSIC  results  from  the  simulation  as  well  as  the  target  setup  dia¬ 
log  window.  Section  5  draws  a  conclusion. 


^Registered  trademark  of  the  Mathworks  Corporation. 
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2.  Extension  of  the  MUSIC  Algorithm  for  Scanning  Antenna 


Superresolution  algorithms,  such  as  the  MUSIC  algorithm,  have  a  special 
structure  that  uses  the  phase  augmentation  in  elements  of  a  linearly  spaced 
sensor  array.  However,  for  a  scanning  antenna,  the  phase  is  almost  constant 
within  the  main  lobe  of  the  transmitted  antenna  beamwidth,  while  the  am¬ 
plitude  of  the  response  varies.  Therefore,  the  returned  signal  does  not  have 
the  same  structure  as  in  the  linearly  or  circularly  spaced  sensor  array. 

The  algorithms  that  use  the  signal  subspace  method  fail  to  resolve  closely 
spaced  targets  within  the  main  lobe  of  the  antenna  beamwidth  because  of 
rank  deficiency.  To  remedy  this  and  to  "decorrelate"  coherent  signals,  Evans 
et  al.  [5]  introduced  the  spatial  smoothing  method,  and  its  proof  was  shown 
in  Shan  et  al.  [6]  and  Friedlander  et  al.  [7].  This  method  was  originally  de¬ 
rived  for  the  linearly  spaced  sensor  array.  We  adapted  the  spatial  smooth¬ 
ing  method  for  a  scanning  radar  antenna,  called  subvector  averaging 
method,  because  the  returned  signal  was  received  in  a  vector  as  the  an¬ 
tenna  stepped  through  discrete  angles  in  the  angular  region  of  interest.  This 
is  a  crucial  step  in  producing  a  correlation  matrix  of  sufficient-rank.  Com¬ 
putational  verification  that  the  subvector  averaging  method  produces  the 
necessary  effect  of  a  rank  increase  in  the  correlation  matrix  with  a  meas¬ 
ured  antenna  pattern  for  point  targets  was  shown  in  other  documented  re¬ 
search  [1].  In  this  section  we  develop  the  SMUSIC  algorithm  using  the 
subvector  averaging  method. 

To  illustrate  the  difference  between  scanning  antenna  and  a  linear  sensor 
array,  figure  1  shows  the  scanning  radar  antenna  beams.  In  this  figure,  the 
main  lobe  of  the  antenna  is  stepping  from  1  to  N  scan  steps  with  a  single 
corner  reflector  at  the  angle  13.  As  the  antenna  step-scans  at  different  angles 
in  an  angular  region  of  interest,  we  obtain  an  input  data  vector  composed 
of  N  elements. 

Assume  that  there  are  two  complex  scatterer-type  targets  which  consist  of  a 
total  of  M  stationary  planewaves  (complex  scatterers  with  their  magnitudes, 
phases,  and  locations).  A  snapshot  in  this  context  is  defined  as  the  radar 
completes  one  step  scan  from  1  to  N.  The  received  signal  for  a  scanning 
antenna  at  the  kth  snapshot  model  with  M  scatterers  is 


'Xik 

SlZc 

Xk  = 

Xlk 

•  1 

XNk 

where  6=  [dy  dj^]^  is  a  vector  of  angles  that  the  antenna  steps  through; 
I3i,  f32f  ■■■f  Pm  are  the  locations  of  scatterers  of  two  complex  targets;  Si^,  where 
z  =  1,  . . .,  M,  is  the  element  of  M  x  1  complex  amplitude  signal  vector  at 
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Figure  1.  Scanning 
antenna  beams. 


Figure  2.  K„  band 
anterma  pattern  at 
35  GHz. 


the  kth  snapshot;  m*.  is  the  N  xl  Gaussian  noise  vector  at  A:th  snapshot,  each 
of  whose  elements  is  independent;  and  a{6, 13, „)  is  an  N  x  1  steering  vector  or 
array  manifold  vector  of  the  measured  scanning  antenna,  shown  in  an  ex¬ 
ample  of  which  is  figure  2*  in  the  direction  m  =  1,  M. 

Thus, 


a{d,l3J 


a{d2-l3m) 


“  Pm) 


A  compact  vector  notation  for  a  scanning  antenna  is 


x=As+u , 


(2) 


(3) 


*This  sum  pattern  is  12"  lense  of  35  GHz  radar  in  which  the  data  were  digitized  from 
Militech  report  (2/19/96). 
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where  xis  aN  x  K  matrix,  K  is  the  total  number  of  snapshots,  sis  aMx  K 
matrix,  u  is  the  N  x  K  noise  matrix,  and  A  is  the  N  x  M  array  manifold 
matrix.  Thus, 


A  =  [a{d,l3^)\aid,l32)\---\a{d,l3i^)]  , 


(4) 


The  important  property  of  representing  A  as  in  equation  (4)  is  that  the  ele¬ 
ments  of  the  vector  a{d,  jS)  for  the  scanning  antenna  are  not  simply  delayed 
versions  of  one  another  as  in  a  uniformly  linear  spaced  sensor  array 


2.1  Subvector  Averaging  Method 


Since  the  SMUSIC  superresolution  algorithm  that  uses  the  signal  subspace 
method  utilizes  the  noise  eigenvector  subspace  to  compute  a 
pseudospectrum,  it  is  important  to  find  the  noise  eigenvector  subspace  from 
the  data  correlation  matrix.  In  this  section,  I  describe  the  subvector  averag¬ 
ing  method  to  compute  the  correlation  matrix.  Given  input  data  x  as  shown 
in  equation  (3),  with  spatially  independent  noise,  the  correlation  matrix  can 
be  found  as 


R  =  E[xx^]  =  ASA^  +0^1 


(5) 


where  E  is  the  ensemble  average  from  the  snapshots,  S  =  £[ss^]  is  the  signal 
covariance  matrix,  and  H  denotes  Hermitian  transpose.  E[uu^]  =  o^I,  where 
is  the  noise  power,  and  I  is  an  identity  matrix. 

The  correlation  matrix  R  in  equation  (5)  is  the  true  correlation  matrix  if  there 
is  almost  an  infinite  number  of  snapshots.  However,  in  reality,  given  a  set  of 
data  with  a  finite  number  of  snapshots  or  even  only  one  snapshot  {k  =  1), 
we  need  to  estimate  the  correlation  matrix  by  using  the  subvector  averag¬ 
ing  method,  since  the  rank  of  R  =  xx^  is  equal  to  1  for  k  =  l. 

Let  us  partition  an  input  data  vector  with  the  subscript  k  being  dropped 
out,  X  =  [xi,  X2,  ...,  with  N  step  scans  into  overlapping  sub  vectors  of 
size  Q,  with  elements  {1,  ...,  Q}  forming  the  first  subvector,  with  elements 
{2,  . . .,  Q  -I- 1}  forming  the  second  subvector,  and  so  on,  as  follows: 


where  P  =  N-Q  +  lis  the  total  number  of  subvectors,  P  <Q. 

Next,  let  the  instantaneous  subcorrelation  matrix  be  x^P^  Finally, 
define  the  estimated  correlation  matrix  to  be  the  average  of  P 
subcorrelation  matrices  as 


(7) 
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This  method  has  not  been  proven  theoretically  for  a  scanning  antenna;  it 
has  been  proven  only  for  a  uniformly  linear  spaced  sensor  array  Other  docu¬ 
mented  research  has  verified  computationally  that  the  subvector  averaging 
method  works  for  scanning  antenna  data  in  [1].  The  theoretical  proof  is 
being  investigated. 

2.2  SMUSIC  Algorithm 

The  matrix  is  of  the  order  Q  x  Q.  If  there  are  two  target  centroids  from 
M  scatterers,  there  will  be  two  (2)  predominant  eigenvalues  from  .  Let 
{A,£}  be  the  eigendecomposition  of  Rxx^  where  the  eigenvalues 
A  =  { ,  A2 , . . . ,  Aq  }  are  arranged  in  decreasing  order  of  size,  with  correspond¬ 
ing  normalized  eigenvectors  E  =  {ei,e2,...,eQ}.  Let  us  partition  these  into 
signal  and  noise  subspace  components  with  Eg  =  {e-i,e2,...,e^}  having  the 
dimension  Q  x  M  and  having  the  dimension  Q  x 

A  .A  A  A  A  A  A 

(Q  -  M),  so  that  E  =  {Eg,Eyj}  and,  similarly,  Ag  =  {Ai,A2,...,Aj\4}  and 
=  {'^M+i/'^M+z/'-'/'^ql/Sothat  A  =  {As,Ajy}.  These  eigenvector  subspaces 
are  orthogonal  to  each  other.  For  this  paper  we  know  there  are  two  targets, 
so  we  can  split  the  signal  subspace  with  dimension  Q  x  2  and  the  noise 
subspace  with  dimension  Q  x  (Q  -  2).  In  a  general  case,  one  could  use  the 
minimum  description  length  technique  [8, 9]  to  estimate  the  number  of  tar¬ 
gets  before  splitting  the  signal  and  noise  subspaces. 

Let  us  introduce  the  a{jS)  vector  to  be  the  effective  array  manifold  vector  for 
all  allowable  angles  yS  in  the  true  signal  space.  In  this  case,  a{jS)  is  the  two- 
way  digitized  magnitude  pattern  of  the  Ka  band  radar.  The  effective  array 
manifold  vector  was  computed  the  same  way  as  in  the  estimate  of  the  cor¬ 
relation  matrix,  i.e.,  the  subvector  averaging  method  with  the  same  number 
of  subvector  averages  was  used  to  compute  the  array  manifold  vector  to 
avoid  the  shifting  in  angular  space.  Note  that  the  number  of  elements  of 
the  estimated  manifold  vector  must  be  the  same  the  number  of  elements  in 
each  column  of  the  noise  eigenvector  subspace.  Since  the  signal  subspace 
contains  information  about  the  angles  of  arrival  from  each  planewave,  the 
array  manifold  vector  from  those  angles  are  also  orthogonal  to  the  vectors 
in  the  noise  subspace.  The  magnitude  of  the  product  of  the  estimated  mani¬ 
fold  vector  and  the  noise  subspace  is  therefore  a  small  number,  i.e., 

\\a^{fi)E^f  =  e  or 

aiP)^Ej^E^a{l3)  =  e  ,  (8) 

where  e  is  a  small  number. 

Equation  (8)  is  a  null  SMUSIC  spectrum.  If  we  normalize  this  equation,  it 
becomes 


ail3)^Ej^E^a{l3) 

a{l3fa{l3) 


(9) 
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The  inverse  of  equation  (9),  i.e.,  the  magnitude  product  between  a  manifold 
vector  from  all  possible  angles  and  the  noise  subspace  matrix,  called  the 
SMUSIC  pseudospectrum,  is 


PmuH^)  = 


for  all  (3  ■ 


(10) 


The  quantity  of  equation  (10)  is  large  for  yS  such  that  a{l3)s  are  in  the  true 
signal  space.  We  see  that  the  peaks  from  this  spectrum  are  more  intuitively 
appealing.  is  the  profile  of  a  one-dimensional  pseudospectrum.  The 

location  of  the  peaks  of  this  spectrum  are  the  estimated  locations  of  the 
target  centriods.  Equation  (10)  only  estimates  the  locations,  not  their  strength. 
There  is  a  separate  task  to  determine  the  strength  of  the  estimated  signals 
by  applying  the  least-square  solution  technique  [10].  This  paper  does  not 
cover  this  task. 


Equation  (10)  is  the  same  way  as  an  ordinary  MUSIC  calculation.  The 
SMUSIC  algorithm  performs  the  search  for  all  possible  angles  yS  in  the  search¬ 
ing  angular  region  to  find  out  where  the  peaks  of  the  SMUSIC  spectrum 
are.  The  searching  method,  which  goes  through  each  computational  loop  as 
shown  in  figure  3  for  each  searching  angle,  consumes  a  lot  of  computa¬ 
tional  time.  Anew  procedure  or  method  has  been  developed,  taking  advan¬ 
tage  of  the  optimized  MATLAB  software,  for  a  fast  computation  of  the 
SMUSIC  spectrum.  The  technique  involves  the  following: 

1.  Pre-calculate  fl(yS)  at  each  angle  for  all  the  searched  angles. 

2.  Store  them  in  a  file  or  a  computer  memory. 

3.  Create  a  searched  array  manifold  matrix,  A(P)  in  which  the  columns 
of  the  matrix  are  the  array  manifold  vector  of  the  searched  angles. 

I  utilized  the  property  of  the  matrix  product  and  the  diagonal  matrix  ma¬ 
nipulation  to  develop  a  fast  technique  to  calculate  the  SMUSIC  spectrum  as 
shown  below 


Pfast  il3)  = - - , 

diag{A{l3fE^E^A{l3)) 


for  all  yS 


(11) 


Erom  equation  (11),  we  see  that  the  numerator  is  the  vector  of  all  the  diago¬ 
nal  elements  of  the  dot  product  of  the  searched  array  manifold  matrix  and 
its  conjugate  transpose.  Since  we  know  A(yS),  the  only  thing  that  we  need  in 
the  denominator  is  to  compute  the  noise  eigenvector  subspace  . 

In  summary,  the  SMUSIC  process  is  shown  in  figure  3.  Note  that  fl(yS)  in  the 
last  block  of  the  diagram  is  replaced  with  A(yS)  for  the  fast  implementation. 
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Figure  3.  Ordinary 
SMUSIC  process  for 
experimental  data 
and  simulated  data. 


Locate 

peaks  of  SMUSIC  spectrum 


lk(|3)lf 


SVA  =  subvector  averaging  method 
SVD  =  singular  value  decomposition 
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3.  Interface  Program  Description 


The  interface  program  was  written  in  the  MATLAB  language,  version  6  or 
R12.  The  MATLAB  code  for  this  program  is  listed  in  the  appendix.  To  run 
this  program,  first  you  need  to  invoke  the  MATLAB  command  window. 
Then,  assuming  that  the  interface  program  is  in  the  current  directory,  type 
smusicinter  at  the  MATLAB  prompt  sign  »  in  the  MATLAB  command 
window.  The  SMUSIC  interface  window  will  appear  on  the  screen.  Figure  4 
shows  the  interface  program  window.  The  window  consists  of  the  plot  of 
the  SMUSIC  result,  and  the  input  parameter  section.  The  input  parameter 
section  includes  three  parts:  Radar  Characteristics,  Target  Characteristics, 
Processing,  and  Run  SMUSIC  and  other  options. 


3.1  Radar  Characteristics 


Figure  4.  Interface 
program  window. 


Radar  characteristics  consist  of 

•  Signal-to-noise  ratio  (SNR):  the  signal  to  noise  ratio  in  dB. 

•  Range  (RANGE)  from  radar  to  target  center:  the  range  from  the  radar 
to  a  referenced  center  point  of  targets  in  meter. 

•  Bandwidth  (BW)  of  the  system:  the  total  system  bandwidth  in  hertz. 
The  operating  frequency  is  35  GHz. 

•  Number  of  frequency  steps  (NFSTEPS):  the  number  of  frequency  steps 
within  the  system  bandwidth.  These  steps  are  evenly  divided. 

To  enter  the  values  of  the  parameters,  use  the  left  mouse  button,  click  in  the 
box,  then  enter  the  value  in  each  box.  For  example:  SNR  =  30  dB,  RANGE  = 
600  m,  BW  =  2  GHz,  and  NFSTEPS  =  51  (evenly  divided  from  34  GHz  to 
36  GHz). 


SMUSIC 
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3.2  Target  Characteristics 

Figure  5  shows  the  target  setup  diagram. 

Target  characteristics  consist  of 

•  Target  separation  (TGTSEP):  the  target  separation  in  term  of  the  an¬ 
tenna  beam  width,  either  Vi,  V3,  V4,  or  Vg  of  a  beamwidth  apart. 

•  Number  of  point  scatterer  (NPTS)  for  each  target:  the  number  of  scat- 
terers  for  each  target.  In  this  version,  we  assume  there  are  two  com- 
plex-scatterer  targets.  Each  target  has  NPTS  scatterers.  Por  example, 
NPTS  =  10  means  there  are  10  scatterers  for  each  target. 

•  Target  orientations  (TGTORI):  the  vector  of  angles  for  orientations  of 
targets  in  degrees.  These  orientations  show  the  rotations  of  the  tar¬ 
gets  around  their  reference  points  (RANGE,  TGTSEP).  TGTORI  =  [0 
45]  means  no  rotation  for  Target  1,  and  Target  2  rotates  45°  around  its 
reference  center  with  respect  to  the  radar  line-of-sight  axis. 

•  Target  coordinations  (TGTSET):  a  flag  which  indicates  that  the  loca¬ 
tions  of  scatterers  of  each  target  whether  they  are  selected  from  a 
predetermined  fixed  point  target  set  or  from  a  random  scatterer  set 
with  their  downrange,  crossrange  generated  from  a  uniform  random 
generator  within  a  reference  target  frame. 

•  Target  amplitude  (TGTAMPL):  a  flag  which  indicates  whether  the 
scatterers  of  each  target  have  fixed  amplitudes  or  random  amplitudes. 
In  the  fixed-amplitude  case,  the  scatterers  have  the  same  amplitude 
of  one.  In  the  random  case,  the  scatterers  have  random  amplitudes 
generated  from  a  uniform  random  generator  from  5  to  15  dBsm. 

•  Target  offset  (TGTOPPSET):  the  down-range  offset  between  the  two 
reference  centers  as  shown  in  figure  5. 


Figure  5.  Target  setup 
diagram. 
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3.3  Processing 


Processing  consists  of 

•  Method  of  computation  (COMP/MAG):  a  flag  which  indicates 
whether  the  SMUSIC  algorithm  uses  complex  or  magnitude  data 
when  it  calculates  the  SMUSIC  spectrum.  For  the  complex  case,  the 
algorithm  uses  the  complex  data  in  the  SMUSIC  calculation  with  the 
magnitude  antenna  pattern  as  shown  in  figure  2.  For  the  magnitude 
case,  the  algorithm  uses  both  the  magnitude  data  and  the  magnitude 
antenna  pattern  in  the  SMUSIC  calculation. 

•  Number  of  subvector  averages  (SVA):  the  number  of  times  that  we 
perform  subvector  averaging  to  estimate  the  data  correlation  matrix. 
For  example,  SVA  =  3  means  that  we  perform  3  partitions  of  the  in¬ 
put  data,  and  form  three  subcorrelation  matrices  and  then  calculate 
the  average  of  those  subcorrelation  matrices.  Usually,  SVA  =  3  is  suf¬ 
ficient  to  compute  the  correlation  matrix  if  we  want  to  estimate  two 
centroids. 

3.4  Run  SMUSIC  and  Other  Options 

•  Run  SMUSIC:  is  the  execution  button  which  starts  the  SMUSIC  cal¬ 
culation. 

•  Plot  target  setup:  plots  the  target  setup  for  the  simulation  after  the 
SMUSIC  calculation  is  done. 

•  Plot  SMUSIC:  plots  the  SMUSIC  spectrum  for  the  simulation  after 
the  SMUSIC  calculation  is  done. 

•  Plot  data  at  frequency  i,  where  i  =  1,  . . .,  NFSTEPS:  allows  a  user  to 
plot  the  data  at  a  particular  frequency. 

•  Info:  shows  where  the  MATLAB  source  code  of  the  interface  pro¬ 
gram  is  located. 

•  ZOOM:  allows  users  to  span  the  figure  at  a  specific  area.  To  do  this, 
first  click  the  ZOOM  button,  and  then  move  the  mouse  cursor  in  the 
desired  area  of  the  figure  and  click  the  left  mouse  button.  If  you  want 
to  zoom  back  out  to  the  normal  size  figure,  click  the  right  mouse  but¬ 
ton. 

•  REFRESH:  clears  all  input  parameters,  figure,  and  screen  command 
window. 

•  QUIT:  exits  the  interface  program. 
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4.  Simulation 


Since  there  are  12  input  parameters  to  the  program,  we  cannot  set  up  all  the 
possible  combinations.  The  result  of  the  simulation  for  one  set  of  the  input 
parameters  is  given  in  table  1. 

Figure  6  shows  the  result  of  the  SMUSIC  spectrum,  plotted  as  SMUSIC  out¬ 
put  (dB)  versus  Azimuth  angle  (degree),  with  the  specified  input  param¬ 
eters.  In  this  figure,  the  two  vertical  lines  on  the  top  of  the  figure  are  the 
estimated  centroids  of  the  two  targets.  The  shaded  area  is  the  antenna 
beamwidth  at  -3  dB  point.  This  figure  clearly  shows  that  the  SMUSIC  algo¬ 
rithm  resolves  the  centroids  of  the  two  targets  accurately  at  V4  of  the  an¬ 
tenna  beamwidth.  The  peaks  of  the  SMUSIC  spectrum  almost  coincide  with 
the  estimated  target  centroids.  Other  simulation  results  can  be  obtained  by 
changing  the  values  of  the  input  parameters. 

Figure  7  shows  the  target  setup  with  the  corresponding  input  parameters 
that  are  shown  in  table  1.  Each  target  has  10  point  scatterers.  Each  scatterer 
has  its  amplitude  and  phase.  The  larger  the  size  of  the  scatterer,  the  larger 
the  amplitude.  The  dash  box  around  each  target  represents  the  reference 
frame  of  a  target  with  a  dimension  of  6  m  in  the  downrange  direction  and 
4  m  in  the  crossrange  direction.  The  center  of  the  reference  box  is  at  (RANGE, 
TGTSEP).  In  this  case,  RANGE  is  at  600  m  down  range  and  TGTSEP  is  V4  of 
a  beamwidth  apart.  The  two  targets  are  facing  directly  to  the  radar,  i.e.,  = 
a2  =  0.  The  shaded  area  in  the  figure  indicates  the  area  covered  by  the  an¬ 
tenna  beamwidth  at  the  -3  dB  point,  about  23.24  m  for  the  antenna 
beamwidth  of  2.2°. 

Pigure  8  is  the  target  response  at  the  frequency  number  30  (35.16  GHz). 
Note  that  there  are  51  different  frequencies  for  this  particular  simulation. 
These  frequencies  are  divided  equally  within  a  2  GHz  interval  with  the  cen¬ 
ter  frequency  at  35  GHz. 

In  summary,  the  interface  tool  is  extremely  useful  for  analyzing  the  perfor¬ 
mance  of  the  SMUSIC  algorithm  for  different  scenarios  with  complex  scat- 
terer-type  targets.  We  can  use  this  tool  to  plan  an  experiment  before  we  can 
actually  conduct  an  experiment. 


Table  1.  Input 
parameters  for 
SMUSIC  simulation. 


Parameter 

Value 

SNR  (dB) 

30 

RANGE  (m) 

600 

BW  (Hz) 

2E09 

NFSTEPS 

51 

TGTSEP  (BW) 

1/4 

NPTS 

10 

TGTORI  (°) 

0 

TGTSET 

"Random" 

TGTAMPL 

"Random" 

TGTOFFSET  (m) 

0 

COMP/MAG 

"Complex" 

SVA 

3 
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Figure  6.  SMUSIC 
result. 


Figure  7.  Target  setup 
window. 
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Figure  8.  Target 
response  at  the 
frequency  number  30 
(35.16  GHz). 
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5.  Conclusions 


This  report  describes  the  extension  of  the  MUSIC  algorithm  to  the  SMUSIC 
for  complex  scatterer-type  targets.  The  extension  showed  clearly  the  differ¬ 
ence  between  the  MUSIC  algorithm  which  used  a  uniformly  linear  sensor 
array  and  the  SMUSIC  algorithm  which  used  a  single  rotatable  sensor.  Also 
described  is  the  graphical  user  interface  program  which  allows  us  to  ana¬ 
lyze  the  performance  of  the  SMUSIC  algorithm.  The  SMUSIC  algorithm 
can  resolve  two  closely  spaced  simulated  complex  scatterer-type  targets  and 
can  estimate  the  locations  of  those  centroids  with  frequency  diversification 
regardless  of  constructive  or  destructive  interference.  The  targets  are  within 
the  -3  dB  beamwidth  of  the  band  radar  antenna.  This  work  shows  that 
we  can  improve  the  angular  resolution  by  at  least  a  factor  of  4.  This  report 
also  describes  a  fast  implementation  of  the  SMUSIC  algorithm,  which  shows 
great  promise  to  be  implemented  in  a  real-time  system. 
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Acronyms 


SMUSIC 

Scan-Multiple  Signal  Classification 

RANGE 

Range  from  the  center  of  the  radar  to  center  of  corner 
reflectors 

BW 

System  bandwidth 

ARE 

Army  Research  Laboratory 

NESTEPS 

Number  of  Erequency  Steps 

TGTSEP 

Target  Separation  in  beam  width 

MATLAB 

Registered  trademark  of  the  Mathworks  Corporation 

NPTS 

Number  of  point  scatterers 

TGTORI 

Target  orientation 

SNR 

Signal  to  noise  ratio 

TGTSET 

Target  set 

TGTAMPL 

Amplitude  of  ears  scatterer  in  the  target  set 

TGTOPPSET 

the  target  of  set  between  the  first  and  second  target 

COMP/MAG 

Computational  selection,  complex  or  magnitude 
computation 

SVA 

Number  of  subvector  averagings 
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Appendix.  MATLAB  Source  Code  of  the  Interface  Program 


%  Filename:  smusicinter . m 
%  Author:  Canh  Ly 
%  Date:  29  November  2000 

%  SMUSICINTER  brings  up  the  SMUSIC  GUI.  There  are  two  sections  for 
%  the  user  to  input  parameters  for  the  simulation.  Section  1  is  the  Radar 
%  Characteristics,  and  section  2  is  the  Target  Characteristics.  The  radar 
%  characteristics  parameters  include  signal-to-noise  ratio  (SNR)  in  dB, 

%  radar  range  (range)  in  meter,  radar  bandwidth  (BW)  in  Hz,  and  number  of 
%  step  frequencies  (nfsteps) .  The  target  characteristic  parameters 
%  include  target  separation  (tgtsep) ,  number  of  point  scatterer  for  each 
%  set  of  a  target  (npts),  target  orientation  (tgtori) ,  target  set,  either 
%  one  point  target  for  each  target  or  multiple  point  targets  for  each  target 
%  set  (tgtset),  and  target  amplitude  (tgtampl) ,  either  point  targets  have 
%  the  same  magnitude  or  random  magnitude  ranged  from  5  dBsm  to  15  dBsm  in 
%  total  number  of  point  scatterers.  The  SMUSIC  result  will  be  printed  in 
%  the  MATLAB  window. 

%  The  GUI  then  plots  the  SMUSIC  pseudo-spectrum  in  dB  versus  azimuth  search 
%  angle . 


nsteps  =  l;value=l;  %  Setup  default  number  of  steps  and  value  for  frequency  selection 
figtag  =  findobj(0,  'tag',  'SMUSIC'); 
if  isempty ( figtag) 

fig  =  figure ( 'Name' ,  'SMUSIC' ,' Units' ,' pixels' ,' Position' ,  [303  188  658  495],... 

'NumberTitle' ,  'off' ,  'resize' ,  'on' ,  'tag' ,  'SMUSIC' ) ; 

end 

set (gcf ,  'DefaultAxesPosition' ,  [.I  .1  .55  .8]); 

figcolor  =  get (gcf,  'color' ) ; 

hndtexthead . Style  =  'text' ; 
hndtexthead . Unit s  =  'normal'; 
hndtexthead . backgroundcolor  =  figcolor; 
hndtexthead . foregroundcolor  =  'blue'; 

hndtextsubhead . Style  =  'text'; 
hndtextsubhead . Units  =  'normal'; 
hndtextsubhead . backgroundcolor  =  figcolor; 
hndtextsubhead . foregroundcolor  =  'black'; 

%  Radar  Characteristic  section 
uicontrol (gcf,  hndtexthead, . . . 

'position',  [.69  .95  .27  .03],... 

' font size' , 10 , ' font weight ','bold',... 

'string' ,  'Radar  Characteristics' ) ; 

uicontrol (gcf,  hndtextsubhead, . . . 

'position',  [.70  .91  .1  .03],... 

'string',  'SNR  (dB) ' ) ; 
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uicontrol (gcf ,  hndtextsubhead, . . . 
'position',  [.85  .91  .1  .03],... 
'string' ,  'RANGE  (m) ' ) ; 

uicontrol (gcf,  hndtextsubhead, . . . 
'position',  [.70  .82  .1  .035],... 
'string',  'BW  (Hz)'); 

uicontrol (gcf,  hndtextsubhead, . . . 
'position',  [.85  .82  .1  .035],... 
'string',  'NFSTEPS'); 

%  Target  Characteristic  section 
uicontrol (gcf,  hndtexthead, . . . 
'position',  [.69  .725  .28  .0385],.. 
' font size'  ,  10 ,  '  font weight ' , 'bold' ,  . 
'string',  'Target  Characteristics') 

uicontrol (gcf,  hndtextsubhead, . . . 
'position',  [.70  .685  .1  .035],... 
'string',  'TGTSEP' ) ; 

uicontrol (gcf,  hndtextsubhead, . . . 
'position',  [.85  .685  .1  .035],... 
'string' ,  'NPTS' ) ; 

uicontrol  (gcf,  hndtextsubhead,  . . . 
'position',  [.70  .595  .1  .035],... 
'string',  'TGTORI'); 

uicontrol  (gcf,  hndtextsubhead,  . . . 
'position',  [.85  .595  .1  .035],... 
'string',  'TGTSET' ) ; 

uicontrol (gcf,  hndtextsubhead, . . . 
'position',  [.70  .505  .1  .035],... 
'string',  'TGTAMPL' ) ; 

uicontrol (gcf,  hndtextsubhead, . . . 
'position',  [.84  .505  .12  .035],... 
'string',  'TGTOFFSET' ) ; 

%  Processing  section 
uicontrol (gcf,  hndtexthead, . . . 
'position',  [.69  .41  .27  .0375],... 
' font size'  ,  10 , ' font weight ' , 'bold' ,  . 
'string' ,  'Processing' ) ; 

uicontrol (gcf,  hndtextsubhead, . . . 
'position',  [.69  .370  .12  .035],... 
'string',  'COMP/MAG'); 

uicontrol (gcf,  hndtextsubhead, . . . 
'position',  [.85  .370  .1  .035],... 
'string',  'SVA' ) ; 
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%  Editable  buttons 
hndedit . Style  =  'edit' ; 
hndedit .Units  =  'normal'; 
hndedit . Fontsize  =  10; 
hndedit . fontweight  =  'bold' ; 

hndpopup . Style  =  'popupmenu' ; 
hndpopup . Units  =  'normal'; 
hndedit . Fontsize  =  10; 
hndedit . fontweight  =  'bold' ; 

inputl  =  uicontrol (gcf ,  hndedit,  'position',  [.7  .86  .115  .05]); 
calll  =  ['SNR  =  str2num (get (geo,  ''string''));']; 
set (inputl,  'callBack' ,  calll); 

input2  =  uicontrol  (gcf ,  hndedit,  'position',  [.85  .86  .10  .05]); 
call2  =  ['range  =  str2num (get (geo,  '' string' '));'] ; 
set(input2,  'callBack',  call2); 

inputs  =  uicontrol  (gcf ,  hndedit,  'position',  [.7  .77  .115  .05]); 
calls  =  [ 'bw  =  str2num (get (geo,  '' string' '));'] ; 
set (inputs,  'callBack',  callS) ; 

input!  =  uicontrol  (gcf ,  hndedit,  'position',  [.85  .77  .10  .05]); 
call!  =  [ 'nsteps  =  str2num (get (geo,  '' string' '));'] ; 
set (input!,  'callBack',  call!); 

tgtsepstr= [ '1/2  BW|l/5  BW|1/!  BW|l/8  BW' ] ; 

tgtsepinput=  uicontrol (gcf ,  hndpopup,  'position',  [.7  .625  .115  .05], 
'string' ,  tgtsepstr, ' backgroundcolor ' , ' white' , ' enable' , ' on' , . . . 
'Interruptible' , ' on' ,  ... 

'foregroundcolor' , 'black' ) ; 

set (tgtsepinput , ' callback' , ' settgtsep' ) ; 

input6  =  uicontrol  (gcf ,  hndedit,  'position',  [.85  .  6S5  .10  .05]); 
call6  =  [ 'npts  =  str2num (get (geo,  '' string' '));'] ; 
set (input6,  'callBack',  call6) ; 

input7  =  uicontrol  (gcf ,  hndedit,  'position',  [.7  .5!5  .115  .05]); 
call7  =  [ 'tgtori  =  str2num (get (geo,  '' string' '));'] ; 
set(input7,  'callBack',  call7); 

tgtsetstr= [ 'Fixed | Random' ] ; 

tgtsetinput=  uicontrol  (gcf ,  hndpopup,  'position',  [.85  .5!5  .10  .05], 
'string' ,  tgtsetstr, ' backgroundcolor' , ' white' , ' enable' , ' on' , . . . 
'Interruptible' , ' on' ,  ... 

'foregroundcolor' , 'black' ) ; 

set (tgtsetinput , ' callback' , ' settgtset' ) ; 

tgtamplstr= [ 'Fixed  I  Random' ] ; 

tgtamplinput=  uicontrol (gcf,  hndpopup,  'position',  [.7  .!55  .115  .05] 
'string' ,  tgtamplstr, ' backgroundcolor' , ' white' , ' enable' , ' on' ,  .  .  . 
'Interruptible' , ' on' ,  ... 

'foregroundcolor' , 'black' ) ; 

set (tgtamplinput , ' callback' , ' settgtampl' ) ; 

inputlO  =  uicontrol (gcf,  hndedit,  'position',  [.85  .!55  .10  .05]); 
calllO  =  [ 'tgtoffset  =  str2num (get (geo,  ''string''));']; 
set (inputlO,  'callBack',  calllO); 


processstr= [ 'Complex | Magnitude' ] ; 

processinput=  uicontrol (gcf,  hndpopup,  'position',  [.7  .325  .115  .05] 
'string' ,  processstr, ' back.gr oundcolor' , ' white' , ' enable' , ' on' , . . . 
'Interruptible' , ' on' ,  ... 

'foregroundcolor' , 'black' ) ; 

set (processinput, ' callback' , ' setprocess' ) ; 

inputl2  =  uicontrol (gcf,  hndedit,  'position',  [.85  .325  .10  .05]); 
calll2  =  [ 'sva  =  str2num (get (geo,  '' string' '));'] ; 
set(inputl2,  'callBack' ,  calll2); 

%  Push  bottons 

hndpush . Style  =  'push' ; 

hndpush . Unit s  =  'normal'; 

calcstart  =  uicontrol (gcf,  hndpush,  'position',  [.76  .26  .13  .05],... 
'string' ,  'Run  SMUSIC' , ' backgroundcolor ' , 'blue' ,  .  .  . 

' foregroundcolor' , ' yellow' , ' callback' ,  'msm35GHzd' ) ; 

plottgtsetup  =  uicontrol  (gcf ,  hndpush,  'position',  [.685  0.20  .15  .05 
'string',  'Plot  target  setup' ,' backgroundcolor' ,' blue' ,.. . 
'foregroundcolor' , 'yellow' , 'callback' ,  'plottgtsetupinter' ) ; 

plotsmusic  =  uicontrol  (gcf ,  hndpush,  'position',  [.845  0.20  .13  .05], 
'string' ,  'Plot  SMUSIC' , 'backgroundcolor' , 'blue' , . . . 

'foregroundcolor' , 'yellow' , 'callback' ,  'plot smuresult ' ) ; 


labelstr= [ ] ; 
if  nsteps  ==  1 

labelstr= [ 'Plot  data 
elseif  nsteps>=2 

temp  =  [  'Plot  data  @ 
for  i  =  2 : nsteps 

tempadd  =  [  'Plot 
labelstr  =  [temp 
temp  =  labelstr; 

end 


@  freq.  1  I ' ] ; 

freq.  1|']; 

data  @  freq. 
tempadd] ; 


num2str (i)  '  I ' ] ; 


end 

labelstr= [ labelstr  'Plot  the  average']; 

plotdata  =  uicontrol (gcf ,  hndpopup,  'position',  [.74  .11  .19 
'string' ,  labelstr, ' backgroundcolor' , ' blue' , ' enable' , ' on' , . . 
'Interruptible' , ' on' ,  ... 

'foregroundcolor' , ' yellow' ) ; 
column  =  get (geo, ' value' ) ; 

set (plotdata, ' callback' ,  'plotdatainter ' ) ; 


.08] ,  .  .  . 


%  Push  bottons 


uicontrol 

(gcf. 

hndpush. 

'position' ,  [ . 

775  . 

085 

.  10 

.05]  , 

'string' , 

'  INFO' ,' callback' ,  'infofile 

' ) ; 

uicontrol 

(gcf. 

hndpush. 

'position' ,  [ . 

685  . 

025 

.  08 

.05]  , 

'string' , 

'ZOOM' , ' callback' ,  'zoom' ) ; 

uicontrol 

(gcf. 

hndpush. 

'position' ,  [ . 

775  . 

025 

.  10 

.05]  , 

'string' , 

'REFRESH' , ' callback' ,  'smusicinit 

' )  ; 

uicontrol 

(gcf. 

hndpush. 

'position' ,  [ . 

885  . 

025 

.  08 

.05]  , 

'string',  'QUIT' ,' callback' ,  'close  all' ) ; 

grid 

clc; 
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%function  settgtsep 
%  Filename:  settgtsep. m 
%  Author  :  Canh  Ly 
%  Date  :  16  January  2001 

sep_value  =  get (geo, ' value' ) ; 
if  sep_value==l 

sep  =  2;  %  For  1/2  BW  apart 
elseif  sep_value==2 

sep  =  3;  %  for  1/3  BW  apart 
elseif  sep_value==3 

sep  =4;  %  for  1/4  BW  apart 
elseif  sep_value==4 

sep  =  8;  %  for  1/8  BW  apart 

end 


%function  settgtset (column) 
%  Filename:  settgtset. m 
%  Author  :  Canh  Ly 
%  Date  :  16  January  2001 

tgtset  =  get (geo, ' value' ) ; 


%function  settgtampl (column) 
%  Filename:  settgtampl. m 
%  Author  :  Canh  Ly 
%  Date  :  16  January  2001 

tgtampl  =  get (geo, 'value'); 


%function  setprocess 
%  Filename:  setprocess. m 
%  Author  :  Canh  Ly 
%  Date  :  16  January  2001 


comp_mag  =  get (geo, ' value' ) ; 
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% function  [ ant , ds , smuout ] =msm35GHzd ( SNR, range, bw, nsteps , sep, npt s , tgtori,tgtset, 
%  tgtampl , comp_mag, sva) 

%  Filename:  msm35GHzd.m  (monopulse  Scan-MUSIC  for  35  GHz) 

%  Author  :  Canh  Ly 
%  Date  :  20  November  2000 

%  Run  SMUSIC  using  monopulse  antenna  for  complex  scatterer  type 
%  targets  for  averaging  all  frequencies. 

%  Syntax:  [ant , ds , smuout ] =msm35GHzd 

%  (SNR,  range,  bw, nsteps,  sep, npts, tgtori, tgtset, tgtampl, comp_mag, sva) ; 

%  where: 


% 

INPUT : 

% 

SNR  : 

power  of  signal  noise  ratio 

% 

bw : 

bandwidth  of  a  system 

% 

nsteps : 

number  of  steps  frequency  from  lower  freq 

to 

upper  freq 

% 

range : 

the 

radial  range  from  radar  to  a  centroid 

of 

complex  target 

% 

sep : 

target  separation  with  respect  to  the  'reference'  centroid 

% 

[2] 

for  1/2  beamwidth  apart 

% 

[4] 

for  1/4  beamwidth  apart 

% 

npts  : 

number  of  scatterer  point  targets 

% 

tgtori : 

target  orientation  [deg] 

% 

tgtset : 

[1] 

for  two  point  targets 

% 

[2] 

for  two  set  of  point  scatterers 

% 

tgtampl : 

[1] 

for  equal  magnitudes  1 

% 

[2] 

for  targets  within  5  dBsm  to  15  dBsm 

% 

comp_mag : 

[1] 

for  complex  process 

% 

[2] 

for  magnitude  process 

% 

sva : 

number  of  subvector  averagings 

OUTPUT : 


% 

ds :  data  set  structure.  It 

consists  of 

following  parameters. 

% 

dsl 

[  10x57 

double ] 

% 

Scatterer 

setup 

for  data  set  #1 

% 

ds2 

[  10x57 

double ] 

% 

Scatterer 

setup 

for  data  set  #2 

% 

ds3 

[  10x57 

double ] 

% 

Scatterer 

setup 

for  data  set  #2  + 

1/16  lambda 

% 

ds4 

[  10x57 

double ] 

% 

Scatterer 

setup 

for  data  set  #2  + 

1/8  lambda 

% 

ds5 

[  10x57 

double ] 

% 

Scatterer 

setup 

for  data  set  #2  + 

3/16  lambda 

% 

ds  6 

[  10x57 

double ] 

% 

Scatterer 

setup 

for  data  set  #2  + 

1/4  lambda 

% 

xradarl 

[  10x1 

double ] 

% 

Scatterer 

radar 

range  (down  range) 

values  of 

% 

data  set 

#1 

% 

xradar2 

[  10x1 

double ] 

% 

Scatterer 

radar 

range  (down  range) 

values  of 

% 

data  set 

#2 

% 

yi 

[  10x1 

double ] 

% 

Scatterer 

cross- 

range  coordinates 

for  data 

% 

set  #1 

% 

y2 

[  10x1 

double ] 

% 

Scatterer 

cross- 

range  coordinates 

for  data 

% 

set  #2 

% 

tgtloc 

[  1x20 

double ] 

% 

Cross-range  (azl 

and  az2)  locations  (deg) 

% 

tgtmag 

[  1x20 

double ] 

% 

Magnitude 

of  scatterers  from  5  to 

15  dBsm 

% 

scantheta 

[209x1 

double ] 

% 

Antenna  step  scan  angle 

% 

Asum 

[209x20 

double ] 

% 

Response 

matrix 

coresponding  to  the  ds. tgtloc 

% 

rcsctgtl 

[600.0722  6.2546 

0.5972]  % 

True  res  centroid  of  target  set  #1, 

% 

[x,y,  thetal] 

% 

rcsctgt2 

[599.2045  -5.6135 

-0.5367] 

%  True 

res  centroid  of  target  set  #2 

% 

[x,y,  theta2] 

% 

SNR 

30  %  Input  signal 

-to-noise 

ratio 

% 

xKsumO 

[209x51 

double ] 

% 

xKsumll6 

[209x51 

double ] 

% 

xKsuml 8 

[209x51 

double ] 

% 

xKsum316 

[209x51 

double ] 
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%  xKsuml4:  [209x51  double] 

%  smuout :  outputs  of  SMUSIC  calculation  (N  x  6) ,  where  N  is  the  length  of  the 
%  search  angle,  and  6  columns  of  the  output.  The  first  column  is  the 

%  searching  angle  vector.  The  second  is  the  SMUSIC  output  for  0  deg. 

%  phase  shift.  The  third  is  the  SMUSIC  output  for  45  deg.  phase  shift. 

%  The  fourth  is  the  SMUSIC  output  for  90  deg.  phase  shift.  The  fifth 

%  is  the  SMUSIC  output  for  135  deg.  phase  shift.  The  sixth  is  the  SMUSIC 

%  output  for  180  deg.  phase  shift. 

%  Example : 

%  [ ant , ds , smuout ] =msm35GHzd ( 30 , 600 , 2E9 , 51 , 4 , 10 , 0 , 2 , 2 , 3 , 0 ) ;  %two  rand.  comp.  sea.  sets 


format  compact 

%  Reset  all  frequencies  according  to  number  of  frequency  steps 


labelstr  =  [ ] ; 

hndpopup . Style  =  'popupmenu' ; 
hndpopup . Units  =  'normal'; 
if  nsteps  ==  1 

labelstr= [ 'Plot  data  @  freq.  1  I ' ] ; 
elseif  nsteps>=2 

temp  =  ['Plot  data  @  freq.  1  I  '  ]  ; 
for  i  =  2 : nsteps 

tempadd  =  [  'Plot  data  @  freq. 
labelstr  =  [temp  tempadd] ; 
temp  =  labelstr; 

end 


num2str (i)  ' I ' ] ; 


end 

labelstr= [ labelstr  'Plot  the  average']; 

plotdata  =  uicontrol  (gef ,  hndpopup,  'position',  [.74  .11  .19 
'string' ,  labelstr, ' backgroundcolor ' , 'blue' , 'enable' , 'on' , . . 
'Interruptible' , ' on' ,  ... 

'foregroundcolor' , 'yellow' ) ; 
column  =  get (geo, ' value' ) ; 

set (plotdata, ' callback' ,  'plotdatainter ' ) ; 


.08] ,  .  .  . 


%load  sumAmat35 .mat  %  gen_sumAmat35 .m  generates  sumAmat35 .mat  with  sumpat  variable 
%seaangle=[-3.0:0.5:-1.2  -1.1:0.005:1.1  1 . 2  :  0 . 5  :  3 . 0 ] ; 

%load  sumAmat 35b . mat  %  gen_sumAmat35 . m  generates  sumAmat35 . mat  with  sumpat  variable 
%seaangle=[-3.0:0.1:-1.2  -1.1:0.05:1.1  1 . 2  :  0 . 1 :  3 . 0 ] ; 


%  Using  array  manifold  vector  within  a  program,  without  loading  sumAmat 3 5c . mat 
load  sumAmat 35c . mat 

seaangle=[-3.0:0.1:-1.2  -1.1:0.02:1.1  1 . 2  :  0 . 1 :  3 . 0 ] ; 

%load  sumAmat35d.mat 

%seaangle=[-3.0:0.1:-1.2  -1.1:0.01:1.1  1.2:0. 1:3.0] 


tic 


[ant,  ds] =gdatad (SNR,  bw,  nsteps,  range, sep, npts, tgtori,tgtset, tgtampl, tgt offset) ; 
if  comp_mag==l 

[RxhatO ] =estrxhatf s (ds . xKsumO , sva) ; 
elseif  comp_mag==2 

[RxhatO ] =estrxhatf s (abs (ds . xKsumO ) , sva) ; 
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end 


%  Compute  eigenvalues  and  eigenvectors 
[eigvalO ,  eigvecO ] =eigvv (Rxhat 0 ) ; 

ntgts  =  2; 

%  Compute  signal  and  noise  subspaces 
[signalO,noiseO]=snspaces (eigvalO , eigvecO , ntgts ) ; 

[nrownoisesum, ncolsnoisesum] =size(noiseO) ; 

%  Fastest  method,  method  #1,  1/30/01 
esum0=noise0 ; 
ncolA=size ( sumpat , 2 ) ; 
sumsuba_sum=0 . 0 ; 
for  i=l : sva 

sumsuba_sum  =  sumsuba_sum+sumpat ( i : nrownoisesum+i-1 ,  :  ) ; 

end 

suba_sum  =  (1/sva) *sumsuba_sum; 

%  suba_sum=sumpat ( 1 : nrownoisesum,  : ) ; 
num  =  diag (suba_sum' *suba_sum) ; 

dem  =  diag (suba_sum' * (esum0*esum0 ' ) *suba_sum) ; 
pmusics0=real (10*logl0 (num . / dem) ) ; 

%  End  of  the  fastest  method  block 

%  Better  method,  method  #2,  1/30/01 
%esum0=noise0 ; 

%ncolA=size (sumpat , 2 ) ; 

%for  i=l:ncolA 

%  Normalized  peak  spectral  MUSIC  in  matrix  computation 
%  suba_sum  =  sumpat ( 1 : nrownoisesum, i ) ; 

%  pmusicsO (i) =real (10*logl0 ( (suba_sum' *suba_sum) / (suba_sum' * (esum0*esum0 ' ) *suba_sum) ) ) 
%  i; 

%end 

%  End  of  the  better  method  block 


%  Slow  method,  method  #3,  1/30/01 
%esum0=noise0 ; 


^sum=0 . 0 ; 

^for  i=l : length (seaangle) 

^%%%%%%%  DO  NOT  DELETE  THIS  BLOCK  %%%%%%%%%%%%% 

‘o  Normalized  peak  spectral  MUSIC  in  vector  computation 
‘o  for  k=l :  ncolsnoisesum 

□  [ scantheta, V] =mono_amv35 (seaangle (i) ) 

□  subv  =  V ( 1 : nrownoisesum) ; 

h  sum  =  sum+abs  (  (esumO  (  :  ,  k)  .  '  )  *subv) ''2; 

h  numsum  =  subv' *subv; 


end 

pmusicO (i) =10*logl0 (numsum/sum) ; 
sum  =  0.0; 
i 


%end 

%  End  of  the  slowest  method  block 


toe 

smuout . seaangle  =  seaangle.'; 
smuout .pmusicsO  =  pmusicsO; 
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%f igure ( 2 ) 

plot (smuout . seaangle, smuout .pmusicsO, ' k' , ' linewidth' ,1.5) 
axis  (  [-3  3  0  max ( smuout . pmusicsO ) +7 . 5 ] ) 

if  bw<  IE 9 

bwstr  =  [ 'BW  =  '  num2str (bw/ 1E6 )  '  MHz,  '  ] ; 

elseif  bw  >=  1E9 

bwstr  =  [ 'BW  =  '  num2str (bw/ 1E9 )  '  GHz,  '  ] ; 

end 

if  tgtset==l 

titlestr= [ 'Two  fixed  scatterer  point  targets']; 

else 

titlestr= [ 'SMUSIC,  2  random  data  sets','  f  =  35  GHz,  ', bwstr,  ...\ 

'  SNR  =  ' , num2str (SNR) , '  dB' ] ; 

end 

title  (titlestr, ' font size' , 8, ' font weight '  ,  'bold'  ) 
xlabel ( 'Azimuth  angle  (\circ)') 
ylabel ( 'dB' ) 
hold  on 

%  Plot  the  'actual'  res  centroid  of  multiple  scatterers 

xlinel= [ds . reset gtl rot ( 3 )  ds . reset gtl rot ( 3 ) ] ; 

ylinel= [max ( smuout . pmusicsO ) +2 . 5  max ( smuout . pmusicsO ) +7.5] ; 

xline2= [ds . rcsctgt2rot  ( 3 )  ds . rcsctgt2rot  ( 3 ) ] ; 

yline2= [max ( smuout . pmusicsO ) +2 . 5  max ( smuout . pmusicsO ) +7.5] ; 

plot  (xlinel, ylinel , ' g-' , ' linewidth' ,1.5) 

plot (xline2, yline2 , ' g-' , ' linewidth' ,1.5) ; 

%  Shade  the  antenna  beamwidth 

xantfill= [-ant .beamwidth/2  -ant . beamwidth/2  ant . beamwidth/2  ant . beamwidth/ 2 ] ; 
yantfill=[0  max ( smuout . pmusicsO ) +7 . 5  max ( smuout . pmusicsO ) +7 . 5  0]; 
tcolor=[0.7  0.7  0.7];p  =  f ill (xantf ill, yantf ill, tcolor) ; 
set (p, 'EaceColor' , tcolor, ' EdgeColor ' , 'None' , ' Eace Alpha' ,0.4); 

hold  off 
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function  [ ant , ds ] =gdatad ( SNR, bw, nsteps , range, sep, npts , tgtori , tgtset , tgtampl, tgtof f set ) 
%  Filename:  gdatac.m  (Generate  input  data  using  35  GHz  antenna  pattern) 

%  Author  :  Canh  Ly 
%  Date  :  20  November  2000 

%  Purpose:  Generates  complex  scatterer  type  targets  with  different 
%  distances  of  the  data  set  #2  without  snapshot  K  input  parameter. 


%  Syntax: 

%  [ ant , ds ] =gdatad ( SNR, bw, nsteps , range, sep, npts , tgtori, tgtset, tgtampl , tgt offset) ; 


%  where: 


% 

bw : 

bandwidth  of  a  system 

% 

nsteps : 

number  of  steps  frequency  from  lower  freq 

to 

upper  freq 

% 

range : 

the 

radial  range  from  radar  to  a  centroid 

of 

complex  target 

% 

sep : 

target  separation  with  respect  to  the  'reference'  centroid 

% 

[2] 

for  1/2  beamwidth  apart 

% 

[4] 

for  1/4  beamwidth  apart 

% 

npts  : 

number  of  scatterer  point  targets 

% 

tgtori : 

target  orientation  [deg] 

% 

tgtset : 

[1] 

for  two  point  targets 

% 

[2] 

for  two  set  of  point  scatterers 

% 

tgtampl : 

[1] 

for  equal  magnitudes  1 

% 

[2] 

for  targets  within  5  dBsm  to  15  dBsm 

%  Example : 

%  [ ant , ds ] =gdatad ( 30 , 2E9, 51 , 600 , 2 , 10 , 45, 2 , 5 ) ;  %  For  35  GHz  radar 


[ant , ds ] =tgt setupd (bw, nsteps , range, sep, npt s , tgtori , tgtset , tgtampl, tgtof f set ) ; 

sigmasq=sigmacon (SNR) ; 

sigman=l/ ( sqrt ( sigmasq) ) ; sigmas=l ; 

i=l; 

for  indexdata=7 :size(ds.dslrot,2) 

tgtphaO  =  dtor ( [ds . dslrot ( : , indexdata) '  ds . ds2rot ( : , indexdata) ' ] ) ;  %  move  or  +0m 
signalO  =  ds . tgtmag . *exp ( j *tgtpha0 ) ; 
noise  =  sqrt ( 1/2 ) *sigman*  ... 

(randn ( size (ds . Asumrot ,  1 ), 1 ) +j *randn ( size (ds . Asumrot , 1 ), 1 )) ;  %for  one  snap¬ 
shot 

xKsumO ( : , i) =ds .Asumrot*signal0 . '  +  noise; 
i  =  i  + 1  ; 

end 

ds.SNR  =SNR; 

ds . xKsum0=xKsum0 ;  %  Target  response  of  superposition  of  two  target  data  sets 
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function  Rxhatf ss=estrxhatf s (x, subband) 

%  Filename:  estrxhat f . bs . m  (Generate  estimated  correlation  matrix  with  forward) 
%  Author  :  Canh  Ly 
%  Date  :  21  March  1999 

%  Generate  estimated  correlation  matrix  with 
%  forward  spatial  smoothing  for  scanning  antenna. 

%  Syntax:  Rxhatfss=estrxhatfs (x, subband)  ; 

%  where 

%  x:  input  data  matrix 

%  Rxhatf ss:  estimated  correlation  matrix  with  forward 

%  spatial  smoothing. 


P=subband; 

[N,  K]  =size  (x)  ; 

M=N-P+1;  %  Number  of  elements  in  each  subarray 
Rxtemp=zeros (M,M) ; 
f lag=0 ; 
for  p=l:P 

Rxtemp=Rxtemp+x (p : M, : ) *x (p :M, : ) ' ; 

M=M+1; 

end 

Rxhatf ss=Rxtemp/ (P*K) ; 


function  [eigval, eigvec] =eigvv (Rxhat ) 


Filename 

Author 

Date 


e  1  gw  .  m 
Canh  Ly 
6/15/95 


%  Purpose 
%  Usage 


Computes  eigenvalues  and  eigenvectors 
[eigval,  eigvec] =eigvv (Rxhat) ; 


[V, D]  =eig (Rxhat ) ; 
lambdai=diag (D) ; 

[eigval, eigvec]  =  sorteig (lambda! , V) ; 
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function  [x,y]  =  sorteig (xin, yin) 

%  Filename:  sorteig. m 
%  Name  :  Canh  Ly 

%  Date  :  6/14/94 

%  Purpose:  Sort  two-column  data.  The  sorting 
%  is  based  on  the  order  of  first  the  column. 

%  Usage  :  [x,y]  =  sorteig (xin, yin) 

%  Modification: 

%  8/29/97:  Moved  function  [x,y]=...  to  the  beginning 

%  of  the  file. 

%  8/29/97:  Initialized  y (...)  =  [ 0 . 0 ] ; 

[Xin, I]  =  sort (xin) ; 

1  =  1'; 

leni  =  length(I); 

y(l:size(yin,l),l:size(yin,2))=[0.0]; 

for  j  =  1 : lenI 
index  =  I(j); 
y ( : ,  j )  =  yin ( : , index)  ; 
end 

X  =  flipud(Xin); 
y  =  fliplr  (y) ; 
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function  [ signal , noise] =sn spaces (eigval , eigvec, nsignals ) 
%  Filename:  snspaces.m 
%  Author  :  Canh  Ly 
%  Date  :  March  6,  1997 


%  Purpose  :  Finds  the  number  of  signal  eigenvalues  and 
%  theis  indeces  and  splits  noise  and  signal 

%  space . 

%  Usage  :  [ signal, noise ] =snspaces (eigval , eigvec, nsignals ) ; 

eigvalm  =  eigval ( 1 : nsignals ) ; 
eigvalindex=l : nsignals ; 


%  Split  signal  and  noise  space  vectors 


ncoleigvec=size (eigvec, 2) ; 
noiseindex=l ; 
signalindex=l ; 

leneigvalindex=length (eigvalindex) ; 
for  j=l : ncoleigvec 

if(j  ==  eigvalindex ( signalindex) ) 
o Ids ignalindex=signal index; 
signal ( : ,  signalindex) =eigvec ( : , j ) ; 
signalindex=signalindex+l ; 
if  ( signalindex>leneigvalindex) 
s ignalindex=o Ids ignal index; 

end 

else 

noise ( : , noiseindex) =eigvec  (  : , j ) ; 
noiseindex=noiseindex+l ; 
end 

end 
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function  [ ant , ds ] =tgtsetupd (bw, nsteps , range, sep, npt s , tgtori,tgtset, tgtampl , tgt offset) 
%  Filename:  tgtsetupd.m 
%  Author  :  Canh  Ly 
%  Date  :  17  November  2000 

%  Purpose:  Generate  random  (x,y)s  scatterer  points  at  a  range  using  34  GHz  antenna 
%  pattern  from  35  GHz  radar. 

%  35  GHz  — >  lambda  =  8.571  mm 

%  dia=12"  or  0.3048  m 

%  spotsize  theta  @  -3  dB=  1 . 22*lambda/dia  =  0.03430  rad  =  34.30  mrad. 

%  spotsize  at  600  m  =  spotsize  theta  *  600  =  20.58  m 

%  Syntax :  [ant ,  ds ] =tgt setupd (bw,  nsteps ,  range,  sep, npt s , tgtori,tgtset, tgtampl , tgt off set) 
%  where: 


% 

bw : 

bandwidth  of  a  system 

% 

nsteps : 

number  of  steps  frequency  from  lower  freq 

to 

upper  freq 

% 

range : 

the 

radial  range  from  radar  to  a  centroid 

of 

complex  target 

% 

sep : 

target  separation  with  respect  to  the  'reference'  centroid 

% 

[2] 

for  1/2  beamwidth  apart 

% 

[4] 

for  1/4  beamwidth  apart 

% 

npts  : 

number  of  scatterer  point  targets 

% 

tgtori : 

target  orientation  [deg] 

% 

tgtset : 

[1] 

for  two  point  targets 

% 

[2] 

for  two  set  of  point  scatterers 

% 

tgtampl : 

[1] 

for  equal  magnitudes  1 

% 

[2] 

for  targets  within  5  dBsm  to  15  dBsm 

%  Example : 

%  [ ant , ds ] =tgt setupd (2.0E9,51,  600,2,  10,  45,2,  1,  5)  ; 

%  Note: 

%  What  would  be  the  strength  of  the  simulated  target  since  it  was  generated  from 
%  5  to  15  dBsm?  1 0*logl 0 ( sum ( (ds . tgtmag) . ^2 ) )  in  dBsm.  "ds.tgtmag"  was  in  voltage. 


format  compact 

%  Calculate  the  frequency  steps,  lambdas  and  lambda  at  the  center  frequency 
freq=35E9; 

speed_of_light  =  299792458; 

freqstep=l in space ( f req- (bw/ 2) , f req+ (bw/ 2) , nsteps ) ; 
lambda=speed_of_light . / freqstep; 
lambdacenter=speed_o f_light /f req; 

tgtori=tgtori+180; 

%  Calculate  the  2D  rotation  matrix 
if  length (tgtori ) ==2 

rotate2dl=[cos (dtor(tgtori(l) ) )  -sin(dtor(tgtori(l) ) ) ;sin(dtor(tgtori (1) ) )  . . . 

cos(dtor(tgtori(l) ) ) ] ; 

rotate2d2=[cos (dtor(tgtori(2) ) )  -sin(dtor(tgtori(2) ) ) ;sin(dtor(tgtori (2) ) )  . . . 

cos(dtor(tgtori(2) ) ) ] ; 
elseif  length (tgtori ) ==1 

rotate2dl= [cos (dtor (tgtori) )  -sin (dtor (tgtori) ) ; sin (dtor (tgtori) )  ... 

cos (dtor (tgtori) ) ] ; 
rotate2d2=rotate2dl ; 

end 
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ds . rotate2dl  =  rotate2dl; 
ds.rotate2d2  =  rotate2d2; 

%  antpat_35GHz . mat  generated  from  antpattern35GHz . m,  with  ant  struct 
load  antpat_35GHz . mat 

spotangle=dtor (ant .beamwidth) ;  %  Convert  to  radian. 

spotrange=range*spotangle;  %  Calculate  spotwidth  on  the  ground  at  the  desired  range. 

ds . spotangle  =  spotangle; 
ds . range  =  range; 

ds.spotrange  =  spotrange; 

%  Setup  the  reference  center  for  the  reference  box 

xrefcenterl  =  range; 

yrefcenterl  =  spotrange/ ( sep*2 ) ; 

xrefcenter2  =  xrefcenterl+tgtoff set; 

yrefcenter2  =  -spotrange/ ( sep*2 ) ; 

ds . xrefcenterl  =  xrefcenterl; 
ds . yrefcenterl  =  yrefcenterl; 
ds . xref center2  =  xrefcenter2; 
ds . yref center2  =  yrefcenter2; 

%  Scatterers  are  set  to  the  dimension  of  a  tank  with  front  view  toward  radar 

ylupper  =  yref center 1+2 . 0 ; 

yllower  =  yref center 1-2 . 0 ; 

y2upper  =  yref center2+2 . 0 ; 

y21ower  =  yref center2-2 . 0 ; 

xupper  =3.0;  %  in  m 

xlower  =  -3.0;  %  in  m 

ds. ylupper  =  ylupper; 
ds. yllower  =  yllower; 
ds.y2upper  =  y2upper; 
ds.y21ower  =  y21ower; 
ds. xupper  =  xupper;  %  in  m 
ds. xlower  =  xlower;  %  in  m 

if  tgtset==l%  Use  the  2  fixed  centroid  points 

disp ( 'Two  scatterer  point  targets  regardless  number  of  targets  from  an  user  input'); 
xl  =  0;  x2  =  0; 

yl  =  yrefcenterl;  y2=yref center2 ; 

%  xl  =  xlower  +  (xupper-xlower )  .*  rand(l,l); 

%  x2  =  xlower  +  (xupper-xlower)  .*  rand(l,l); 

%  yl  =  yllower  +  (ylupper-yllower )  .*  rand(l,l); 

%  y2  =  y21ower  +  (y2upper-y21ower )  .*  rand(l,l); 

titlestr= [ 'Two  fixed  scatterer  point  targets']; 
elseif  tgtset==2 

%  Generate  2  random  data  sets  within  reference  boxes 
xl  =  xlower  +  (xupper-xlower)  .*  rand (npts , 1 ) ; 

x2  =  xlower  +  (xupper-xlower)  .*  rand (npts , 1 ) ; 

yl  =  yllower  +  (ylupper-yllower)  .*  rand (npts , 1 ) ; 

y2  =  y21ower  +  (y2upper-y21ower )  .*  rand (npts , 1 ) ; 

titlestr= [ 'Two  random  scatterer  sets  for  35  GHz  radar']; 
end 

%  Rotate  the  data  set  around  the  center  of  the  reference  point,  SP :  Set  Point 
SPl  =  [xl'  ;  (yl-yrefcenterl ) ' ] ;  SP2  =  [x2'  ; (  y2-yref center2 ) ' ] ; 

SPlr  =  rotate2dl*SPl;  SP2r  =  rotate2d2*SP2; 


33 


xlr  =  ( SPlr ( 1 ,  : )  ) ' ;  x2r  =  ( SP2r  ( 1 ,  : ) +tgtof f set ) ' ; 

ylr  =  (SPlr  (2,  : ) ) ' +Yrefcenterl;  y2r  =  (SP2r  ( 2 ,  : ) ) ' +yref center2 ; 

[thrlr, rlr]=cart2pol (xlr, ylr-yrefcenterl ) ; 

[thr2r,  r2r]=cart2pol (x2r, y2r-yrefcenter2 ) ; 

%  Rotational  data  set 

xradarlr=rangetxlr; 

xradar2r=range+x2r; 

%xradar2r=xradar2r+ ( ( 4 / 1 6 ) *lambdacenter ) ;  %  Change  phase  for  two  points 
%  Use  only  azl  and  az2  for  the  target  locations  of  two  scatterer  data  sets 
[thetalr,  Rlr] =cart2pol (xradarlr, ylr) ; 
azlr=rtod (thetalr ) ; 

[theta2r,  R2r] =cart2pol (xradar2r, y2r) ; 
az2r=rtod (theta2r ) ; 


%  2-way  phase  shift,  Deltatheta=  (4*pi/lainbda) ’^Delta_range 
for  i=l : length (lambda) 

philr (:,i)=rtod(mod( ( ( (4*pi) /lambda (i) ) *Rlr) ,2*pi) ) ; 
phi2r ( : , i)=rtod(mod( ( ( (4*pi) / lambda (i) ) *R2r)  ,2*pi) )  ; 

end 

if  tgtampl==l 

ampll=ones ( length (xl ), 1 ) ;  %  All  scatterers  have  the  same  strength. 
ampl2=ones ( length (xl ), 1 ) ;  %  All  scatterers  have  the  same  strength, 
elseif  tgtampl==2 


0  Use  this  block  if  the  statistics  toolbox  is  available. 

0  Generate  uniform  random  from  5  dBsm  to  15  dBsm 

%  5  =  10*logl0 (x^2)  ->  x^2=10^ (5/10)=3.16 

%  the  power,  PiSdBsm,  of  each  scatterer 

%  is  3.16/#  of  scaterers 

%  ViSdBsm  =  sqrt (PiSdBsm) ;  each  scatterer 

%  in  voltage,  Vi 

%  15  =  10*logl0  (x''2)  ->  x''2  =  10^  (15/10)  =31 . 6 

%  the  power,  PilSdBsm,  of  each  scatterer 

%  is  31.6/#  of  scaterers 

%  VilSdBsm  =  sqrt (PilSdBsm) ;  each  scatterer 

%  in  voltage,  Vi 

%  For  example:  for  10  point  scatterers 
%  ViSdBsm  =  0.5623;  and  VilSdBsm  =  1.7783. 

^ampll=unif rnd (ViSdBsm, VilSdBsm, length(xl) , 1) ; 

□  Generate  uniform  random  from  5  dBsm  to  15  dBsm 

%  5  =  10*logl0  (x^2)  ->  x^2  =  10'' (5/10)=3.16 

%  the  power,  PiSdBsm,  of  each  scatterer  is 

%  3.16/#  of  scaterers 

%  ViSdBsm  =  sqrt (PiSdBsm) ;  each  scatterer 

%  in  voltage,  Vi 

%  15  =  10*logl0  (x^2)  ->  x^2  =  10^ (15/10) =31 . 6 

%  the  power,  PilSdBsm,  of  each  scatterer 

%  is  31.6/#  of  scaterers 

%  VilSdBsm  =  sqrt (PilSdBsm) ;  each  scatterer 

%  in  voltage,  Vi 

^ampl2=unif rnd (ViSdBsm, VilSdBsm, length (x2) , 1) ; 
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%  Target  power  between  5  to  15  dBsm 
ViSdBsm  =  sqrt ( ( 10 ^ ( 5/ 10 ) ) /npts ) ; 
VilSdBsm  =  sqrt  (  ( 1 0'' ( 15/ 10  )  ) /npts  )  ; 


ampll  =  ViSdBsm  +  (Vil5dBsm-Vi5dBsm)  . rand  ( length  (xl ),  1 )  ; 
ampl2  =  ViSdBsm  +  (VilSdBsm-ViSdBsm)  . rand  ( length  (x2  ),  1 )  ; 

end 

dslr  =  [Rlr  xlr  ylr  rlr  azlr  ampll  philr] ; 
ds2r  =  [R2r  x2r  y2r  r2r  az2r  ampl2  phi2r] ; 


ds . dslrot=dslr; 

%  Scatterer 

setup  for 

data 

set  #1 

ds . ds2rot=ds2r; 

%  Scatterer 

setup  for 

data 

set  #2 

ds . xradar Irot  = 

#1 

ds . xradar2rot  = 

xradarlr;  % 

Scatterer 

radar 

range 

(down 

range, 

600 

m) 

values 

of 

data 

set 

xradar2r;  % 

Scatterer 

radar 

range 

(down 

range, 

600 

m) 

values 

of 

data 

set 

#2 

ds.xlrot  =  xlr;  %  Rotational  scatterers 
ds.x2rot  =  x2r; 
ds.ylrot  =  ylr; 
ds.y2rot  =  y2r; 

%  Cross-range  (azl  and  az2)  locations  in  polar  coordinate  (deg) 
ds.tgtlocrot  =  [ds . dslrot ( : , 5) '  ds . ds2rot ( : , 5 ) ' ] ; 

ds.tgtmag  =  [ds . dslrot (:, 6) '  ds . ds2rot ( : , 6) ' ] ;  %  Magnitude  of  scatterers 

[ scantheta,  Asumr ] =mono_amv35 (ds . tgtlocrot ) ;  %  Step  scan  angle  and  Response  steering  ma¬ 
trix 

ds.scantheta  =  scantheta;  %  Step  scan  angle 

ds.Asumrot  =  Asumr;  %  Response  matrix  coresponding  to  the  ds.tgtloc 

[xclr,yclr,t he talr]=rcscentroid(ds. dslrot ( : , 1) ,ds . dslrot ( : , 3 ) , ds . dslrot ( : ,  6 ) )  ; 
rcsctgtlr= [xclr  yclr  thetalr]; 

[xc2r,yc2r,theta2r]=rcscentroid(ds.ds2rot ( : , 1) ,ds . ds2rot ( : , 3 ) , ds . ds2rot ( : , 6 ) ) ; 
rcsctgt2r= [xc2r  yc2r  theta2r]; 

ds . rcsctgtlrot=rcsctgtlr;  %  Adding  to  a  structure  with  res  centroid  target  1 
ds . rcsctgt2rot=rcsctgt2r;  %  Adding  to  a  structure  with  res  centroid  target  2 
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function  [ scantheta, sumpat ] =inono_amv35 (tgtloc) 

%  Filename:  mono_amv35.m  (Array  Manifold  Vector  for  monopulse  antenna) 
%  Author  :  Canh  Ly 
%  Date  :  18  September  2000 

%  Generates  array  manifold  vectors  for  35  GHz  monopulse  antenna. 

%  Syntax:  [scantheta,  sumpat] =mono_amv35 (tgtloc) ; 

%  Example : 

%  [ scantheta,  sumpat ] =mono_amv35 ( [-0 . 5  0.5]); 


[beamwidth, outsum] =monol2in35GHz;  %  outsum= [xants  yants  yantsvol]; 

index=f ind (outsum ( : , 1 ) >=-4  &  outsum (:, 1 ) <=4 ) ; 
scantheta=outsum (index, 1) ; 
yant=outsum ( : , 3 ) ; 

for  tgtlocindex=l : length (tgtloc) 

for  scanindex=l : length (scantheta) 

offtheta  =  scantheta(scanindex)-tgtloc(tgtlocindex); 
index=locate (outsum ( : , 1 )  ,  offtheta)  ; 
sumpat (scanindex, tgtlocindex)  =  yant  (index) ; 

end 

end 


function  [beamwidth, outsum] =monol2in35GHz 

%  Filename:  monol2in35Ghz .m  (Sum  pattern  of  35  GHz  monopulse  antenna) 
%  Author  :  Canh  Ly 
%  Date  :  15  September  2000 

%  Generates  the  sum  pattern  of  the  12"  lense  of  35  GHz  radar. 

%  Syntax:  [beamwidth, outsum] =monol2in35GHz; 

%  Example : 

%  [beamwidth, outsum] =monol2in35GHz ; 

%  plot (outsum ( : , 1 ) , outsum ( : , 2 ) ) ; grid 


load  antpat35GHz.dat  %antpat35GHz.dat  data  file  was  digitized  from  Militech  (2/19/96) 
antpat=antpat35GHz ( 1 : 3 : size (antpat35GHz , 1 ) , : ) ; 

[xants ,  yants ] =ant smooth (antpat ( : , 1 ) , antpat  (  : , 2 ) , 10 ) ; 
yants=yants-max (yants ) ; 
yantsvol  =  10 . ^ (yants . /2  0 ) ; 

index=find ( (yants) >-3 . 1  &  (yants ) <-2 . 88 ) ; 

outsum= [xants'  yants'  yantsvol']; 
beamwidth= (xants ( index (2 ) ) -xants ( index ( 1 ) ) ) ; 


function  [xants , yants ] =ant smooth (xantus , yantus , pt s_smooth) 
%  Filename:  antsmooth.m  (Smooth  Antenna  Pattern) 

%  Name  :  Canh  Ly 
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%  Date  :  14  September  2000 

%  This  function  smooths  the  antenna  patterns. 

%  Syntax: 

%  [xants , yants ] =antsmooth (xantus , yantus , pt s_smooth) ; 


n=length (xantus ) ; 
index=l ; 

for  k=l : n- (pts_smooth-l ) 

xants ( index) = smooth (xantus (k : k+pt s_smooth-l ) ) ; 
yants ( index) = smooth (yantus (k : k+pt s_smooth-l ) ) ; 
index=index+l ; 

end 


function  [xs ] =smooth (x) 

%  Filename:  smooth. m  (Smoothing  computation) 
%  Name  :  Canh  Ly 

%  Date  :  Jan.  8,  1998 

%  This  function  smooths  a  data  vector. 

%  Syntax:  [xs ] =smooth (x) ; 

n=length  (x) ; 

coef=diag (fliplr (pascal  (n) ) ) ; 
xs  =  sum  (x .  ’^coef )  /  ( sum  (coef )  )  ; 
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function  [xc, yc, theta] =rcscentroid (x, y, ampl) 

%  Filename:  centroid. m 
%  Author  :  Canh  Ly 
%  Date  :  17  October  2000 

%  Purpose:  Computes  the  centroid  of  the  collection  of  scatterers  with 
%  their  amplitudes . 

%  Reference:  R.  Ellis  and  D.  Gulick, "Calculus  with  Analytic  Geometry," 

%  Second  Edition,  Harcourt  Brace  Jovanovich  (HBJ) ,  pp.  444-445,  1982 

%  Test  case  from  R.  Ellis  et .  al . ,  pp .  445 
%  x=[l  -3  -1  0];y=[2  1-23]; 

%  ampl=[2  457];  %  needed  to  edit  this  file  to  change  this  mass  vector,  m. 

%  [xc, yc, theta] =rcscentroid (x, y, ampl) ; 

%  Answer:  xc  =  -5/6  =  -0.8333;  yc=19/18=l . 0556 

%  Syntax:  [xc,  yc,  theta] =rcscentroid (x, y, ampl) ; 

%  x:  X  coordinates  of  the  object 

%  y:  y  coordinates  of  the  object 

%  m:  amplitude  values  of  the  object  at  locations  (x,y),  m  in  linear  not  in  dB 


%  Assume  all  amplitude  of  scatterers  are  the  same  strength. 
%ampl  =  ones  ( size  (x, 1 ), 1 ) ; 

%amp= [2  4  5  7 ] ; 

%  Calculate  the  moments  of  given  scatterers. 

Mx  =  sura (ampl . *x) ; 

My  =  sum (ampl . *y) ; 

%  Calculate  centroid  point 
m  =  sum (ampl) ; 
xc  =  Mx/m; 
yc  =  My/m; 

%  Now  convert  it  to  angle  in  cross  range  direction 
if  xc>=0 

theta=rtod (atan (yc/xc)  )  ; 
elseif  xc<0 

theta=rtod (pi+at an (yc/xc) ) ; 

end 
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%  Filename:  plottgtsetupinter .m 
%  Author  :  Canh  Ly 
%  Date  :  9  January  2001 

%  Purpose:  Plot  target  setup  for  the  SMUSIC  simulation. 


pltindex  =  1 ; 
for  i  =  1 : length (ds . tgtmag) /2 

if  ds . tgtmag ( i ) *1 00  >=  50  &  ds . tgtmag ( i )* 1 00  <=  76 


msizel  =  2; 

elseif  ds . tgtmag ( i ) *1 00  >  76 
msizel  =  3; 

elseif  ds . tgtmag ( i ) *1 00  >  101 
msizel  =  4; 

elseif  ds . tgtmag ( i ) *1 00  >  127 
msizel  =  5; 

elseif  ds . tgtmag ( i ) *1 00  >  152 
msizel  =  6; 

elseif  ds . tgtmag ( i ) *1 00  >  178 
msizel  =  7; 


ds . tgtmag ( i ) *100  <=  101 
ds . tgtmag ( i ) *100  <=  127 
ds . tgtmag ( i ) *100  <=  152 
ds . tgtmag ( i ) *100  <=  178 


end 

plot (ds . xradarlrot ( i ) , ds . ylrot ( i ) , ' r>' , ' MarkerEdgeColor ' , ' r ' ,  . . . \ 

'MarkerFaceColor ' , ' r' ,  . . . \ 

'MarkerSize' , msizel) 
if  pltindex==l ; 
hold  on 
pltindex  =  0; 

end 

if  ds . tgtmag ( i+ (length (ds . tgtmag) /2 )) *1 00  >=  50  & 

ds . tgtmag ( i+ ( length (ds . tgtmag) /2 )) *100  <=  76 
msize2  =  2; 

elseif  ds . tgtmag ( i+ (length (ds . tgtmag) /2 )) *100  >  76  & 

ds . tgtmag ( i+ ( length (ds . tgtmag) /2 )) *100  <=  101 

msize2  =  3; 

elseif  ds . tgtmag ( i+ (length (ds . tgtmag) /2 )) *100  >  101  & 

ds . tgtmag ( i+ (length (ds . tgtmag) /2 )) *100  <=  127 

msize2  =  4; 

elseif  ds . tgtmag ( i+ (length (ds . tgtmag) /2 )) *100  >  127  & 

ds . tgtmag ( i+ (length (ds . tgtmag) /2 )) *100  <=  152 

msize2  =  5; 

elseif  ds . tgtmag ( i+ (length (ds . tgtmag) /2 )) *100  >  152  & 

ds . tgtmag ( i+ ( length (ds . tgtmag) /2 )) *100  <=  178 

msize2  =  6; 

elseif  ds . tgtmag ( i+ (length (ds . tgtmag) /2 )) *100  >  178 
msize2  =  7; 


end 


plot (ds . xradar2rot (i) ,ds.Y2rot (i) , 'r>' , ' MarkerEdgeColor' , ' r ' ,  . . . \ 

'MarkerEaceColor ' , ' r' ,  . . . \ 

'MarkerSize' ,msize2) 

end 

xorg=0 ; Yorg=0 ; 
scanthe=dtor ([-4:0.05:4]); 
for  i=l : length  (scanthe) 

xscan(i)  =  (ds.range+50) *cos (scanthe  (i) ) ; 
yscan(i)=(ds.range+50) *sin(scanthe(i) ) ; 
xline=[xorg  xscan(i)]; 
yline=[yorg  yscan(i)]; 
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end 

lineanglexl= [ 0  (ds . range+50 ) ] ; lineanglex2= [ 0  (ds . range+50 ) ] ; 
spotrange650= (ds . range+50 ) *ds . spotangle; 
lineangleyl= [ 0  spotrange650 /2 ] ; 

lineangley2= [ 0  ~spotrange650/2 ] ;  %the  spotrange  at  650  m  with  94.44  GHz 

%plot (lineanglexl , lineangleyl , ' b' , lineanglex2 , lineangley2 , ' b' , ' linewidth' ,1.5) 

xfill=[0  (ds . range+50 )  (ds . range+50 )  0]; 

yfill=[0  spotrange650 /2  -spotrange650 /2  0]; 

tcolor=[0.7  0.7  0.7];p  =  fill (xfill, yfill, tcolor) ; 

set (p,  'FaceColor'  , tcolor, ' EdgeColor ' , 'None' , ' Face Alpha' ,0.4); 


%shading  interp; 


%  Plot  true  res  centroid  points 

plot (ds.rcsctgtlrot (1) , ds.rcsctgtlrot (2) , 'go'  ,  '  MarkerEdgeColor '  ,  '  g'  ,  . . . \ 

'MarkerFaceColor '  ,  '  g'  ,  .  .  .  \ 

'MarkerSize' , 5) 

plot (ds.rcsctgt2rot (1) , ds.rcsctgt2rot (2) , 'go' ,  '  MarkerEdgeColor'  ,  '  g'  ,  . . . \ 

'MarkerFaceColor'  ,  '  g'  ,  .  .  .  \ 

'MarkerSize' , 5) 


%  Rotate  boundary  box 

xylbl  =  [ds.xlower  ds . xlower; ds . yllower-ds 
xy2bl  =  [ds.xlower  ds . xupper ; ds . y lupper-ds 
xy3bl  =  [ds. xupper  ds . xupper ; ds . ylupper-ds 
xy4bl  =  [ds. xupper  ds . xlower ; ds . yllower-ds 
xylblr  =  ds  .  rotate2dl’'xylbl; 
xy2blr  =  ds  .  rotate2dl’'xy2bl; 
xy3blr  =  ds  .  rotate2dl’'xy3bl; 
xy4blr  =  ds  .  rotate2dl’'xy4bl; 


yrefcenterl  ds 
yrefcenterl  ds 
yrefcenterl  ds 
yrefcenterl  ds 


ylupper-ds .yrefcenterl]  ; 
ylupper-ds .yrefcenterl]  ; 
yllower-ds . yrefcenterl ]  ; 
yllower-ds . yrefcenterl ]  ; 


xylb2  = 
xy2b2  = 
xy3b2  = 
xy4b2  = 
xylb2r  = 
xy2b2r  = 
xy3b2r  = 
xy4b2r  = 


[  ds 
[  ds 
[  ds 
[  ds 
ds 
ds 
ds 
ds 


.xlower  ds . xlower; ds . y21ower-ds 
.xlower  ds . xupper ; ds . y2upper-ds 
.xupper  ds . xupper ; ds . y2upper-ds 
.xupper  ds . xlower; ds . y21ower-ds 
. rotate2d2*xylb2 ; 

.  rotate2d2’^xy2b2 ; 

.  rotate2d2’^xy3b2; 

.  rotate2d2’^xy4b2 ; 


. yref center2 
.  yref center2 
. yref center2 
. yref center2 


ds . y2upper-ds . yref center2 ] ; 
ds . y2upper-ds . yref center2 ] ; 
ds . y21ower-ds .yrefcenter2] ; 
ds . y21ower-ds .yrefcenter2] ; 


plot (xylblr ( 1 , 
plot  (xy2blr ( 1 , 
plot  (xy3blr ( 1 , 
plot  (xy4blr ( 1 , 
plot  (xYlb2r ( 1 , 
plot  (xy2b2r ( 1 , 
plot  (xY3b2r ( 1 , 
plot  (xY4b2r ( 1 , 


) +ds.xrefcenterl, xylblr (2, :)+ds.yrefcenterl,'b:',' linewidth' ,1.5) 

) +ds.xrefcenterl, xy2blr (2, :)+ds.yrefcenterl,'b:',' linewidth' ,1.5) 

) +ds.xrefcenterl, xy3blr (2, :)+ds.yrefcenterl,'b:',' linewidth' ,1.5) 

) +ds.xrefcenterl, xy4blr (2, :)+ds.yrefcenterl,'b:',' linewidth' ,1.5) 

) + (ds . xref center 2 ) , xylb2r (2, :)+ds.yrefcenter2,'b:',' linewidth' ,1.5) 
) + (ds . xref center2 ) , xy2b2r (2, :)+ds.yrefcenter2,'b:',' linewidth' ,1.5) 
) + (ds . xref center 2 ) , xy3b2r (2, :)+ds.yrefcenter2,'b:',' linewidth' ,1.5) 
) + (ds . xref center2 ) , xy4b2r (2, :)+ds.yrefcenter2,'b:',' linewidth' ,1.5) 


hold  off 


if  ds . range==600 . 0 

axis ( [ds . range-20  ds.range+20 
elseif  ds . range==1000 . 0 

axis ( [ds . range-30  ds.range+30 
elseif  ds . range==1500 . 0 

axis ( [ds . range-40  ds.range+40 
elseif  ds . range==2000 . 0 

axis ( [ds . range-50  ds. range+50 


-20 

20 

-30 

30 

-40 

40 

-50 

50 

40 


end 

xlabel ( 'Down  range  (m) ' ) 
ylabel (  'Cross  range  (m) ' ) 

%title  (  'Target  setup'  ,  '  fontsize' , 12, ' fontweight' , 'bold' ) 
axis  square 
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%  Filename:  plotsmuresult . m 
%  Author  :  Canh  Ly 
%  Date  :  12  January  2001 

%  Purpose:  Plot  the  SMUSIC  simulation  result. 

plot ( smuout . seaangle, smuout . pmusicsO , ' linewidth' ,2.0) 
axis  (  [-3  3  0  max ( smuout . pmusicsO ) +7 . 5 ] ) 

if  bw<  1E9 

bwstr  =  [ 'BW  =  '  num2str (bw/ 1E6 )  '  MHz,  ' ] ; 

elseif  bw  >=  1E9 

bwstr  =  [ 'BW  =  '  num2str (bw/ 1E9 )  '  GHz,  '  ] ; 

end 

title  (titlestr, ' font size' , 8, ' font weight '  ,  'bold'  ) 
xlabel ( 'Azimuth  angle  (\circ)') 
ylabel ( 'dB' ) 
hold  on 

%  Plot  the  'actual'  res  centroid  of  multiple  scatterers 

xlinel= [ds . reset gtl rot ( 3 )  ds . reset gtl rot ( 3 ) ] ; 

ylinel= [max ( smuout . pmusicsO ) +2 . 5  max ( smuout . pmusicsO ) +7.5] ; 

xline2= [ds . rcsctgt2rot  ( 3 )  ds . rcsctgt2rot  ( 3 ) ] ; 

yline2= [max ( smuout . pmusicsO ) +2 . 5  max ( smuout . pmusicsO ) +7.5] ; 

plot  (xlinel, ylinel , ' g-' , ' linewidth' ,1.5) 

plot (xline2, yline2 , ' g-' , ' linewidth' ,1.5) ; 

%  Shade  the  antenna  beamwidth 

xant f ill= [ -ant . beamwidth/2  -ant . beamwidth/2  ant . beamwidth/ 2  ant . beamwidth/2 ] ; 

yantfill=[0  max ( smuout . pmusicsO ) +7 . 5  max ( smuout . pmusicsO ) +7 . 5  0]; 

tcolor=[0.7  0.7  0.7];p  =  fill (xantf ill, yantf ill, tcolor) ; 

set (p, 'EaceColor' , tcolor, ' EdgeColor ' , 'None' , ' Eace Alpha' ,0.4); 

hold  off 
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%  Filename:  infofile.m 
%  Author:  Canh  Ly 
%  Date:  17  January  2001 

set (gcf , ' pointer' , ' arrow' ) 

ttlStr  =  get (gcf ,' name' ) ; 
hlpStrl  =  [ . . . 

'  This  graphical  user  interface  tool  lets  you  run  and  analyze  the 
'  SMUSIC  algorithm.  It  allows  you  to  enter  a  combination  of 
'  parameters  necessary  for  the  simulation  to  analyze  the 

'  performance  of  the  algorithm.  There  are  three  different  sections 

\ 

'  1.  Radar  Characteristics: 

'  In  this  sections  you  can  set  Signal-to-noise  ratio  (SNR)  in  dB, 

'  the  range  (RANGE)  a  distance  between  the  center  of  the 
'  target  and  the  radar  in  meter,  the  total  sytem  bandwidth  (BW)  in 
'  Hertz,  and  the  number  of  step  frequencies  (NFSTEPS)  within 

'  the  bandwidth  by  using  the  controls  on  the  right  of  the  figure. 

\ 

'  2.  Target  Characteristics: 

'  In  this  sections  you  can  set  the  target  separation  (TGTSEP)  in 
'  term  of  the  antenna  beamwidth,  the  number  of  point  scatterers  fo 
'  each  target  (NPTS)  the  vector  of  target  orientations  (TGTORI)  in 
'  degrees,  the  target  is  whether  selected  from  a  random 
'  (Random)  scatterer  set  or  a  fixed  (Fixed)  point  target  set 
'  (TGTSET) ,  the  target  in  which  the  scatterers  of  each  target  have 
'  the  same  amplitudes  (Fixed)  or  random  (Random)  amplitudes 
'  (TGTAMPL) ,  and  the  down  range  offset 

'  (TGTOFFSET)  in  meter  from  the  reference  center  of  one  target 
'  to  the  reference  center  of  the  other  target  by  using  the 
'  controls  on  the  right  of  the  figure. 

hlpStr2  =  [ . . . 

'  3.  Processing: 

'  In  this  sections  you  can  select  whether  the  process  using  the 
'  complex  or  magnitude  process  (COMP /MAG)  in  the 
'  SMUSIC  calculation.  Note  that  in  the  complex  process, 

'  the  array  manifold  vector  is  using  the  magnitude  antenna 
'  pattern  and  complex  data.  In  the  magnitude  process,  both  the 
'  array  manifold  vector  and  the  data  are  using  magnitude.  You 
'  also  can  set  the  number  of  subvector  averages  (SVA)  in  the 

'  estimated  correlation  matrix  calculation. 

\ 

'  Filename:  smusicintera .m 
myFig  =  gcf; 

helpsmusic (ttlStr, hlpStrl, hlpStr2) ; 
return 

\end{ verbatim} 

\pagebreak 

\noi 

\baselineskip=Opc 
\lineskip=0pc 
\lineskiplimit=Opc 
\footnotesize 
\begin{ verbatim} 


%  Filename:  smusicinit.m 
%  Author:  Canh  Ly 
%  Date:  17  January  2001 

clc; 

clf; 

clear  all; 
smusicinter; 
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13.  Abstract  (cont'd) 


of  the  simulation  and  to  demonstrate  the  SMUSIC  algorithm  for  the  complex  scatterer-type 
targets,  I  developed  a  graphical  interface  tool  to  run  and  analyze  the  SMUSIC  algorithm.  The 
tool,  which  is  written  in  the  MATLAB  language,  allows  an  end  user  to  analyze  the  performance 
of  the  SMUSIC  algorithm  by  entering  a  combination  of  parameters  necessary  for  the  simulation. 
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